#include "config.h"

#include <clx/libECM.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_blas.h>

#include <math.h>

#include "lstatfit.h"

#define ERR(fmt, ...)	ECMERR("ERR: lstatfit: " fmt, ## __VA_ARGS__)

char	*inputpath;

static char	*output = "output.txt";

extern ECMBOOL reader_jumpdist(char *buf, ecm_double *pxvalue, ecm_double *pyvalue);

typedef enum {
	JUMPDIST_ALL = 1,
	JUMPDIST_FORWARD,
	JUMPDIST_BACKWARD,
	DOFIT_TYPE_MAX = JUMPDIST_BACKWARD
} dofit_type_t;

typedef struct {
	const char	*desc;
	const char	*logname;
	ecm_uint32	xdiff;
	data_reader_t	reader;
} dofit_info_t;

static dofit_info_t	dofit_infos[] = {
	{ "jumpdist_all", "jumpdist_all", 1, reader_jumpdist },
	{ "jumpdist_forward", "jumpdist_forward", 1, reader_jumpdist },
	{ "jumpdist_backward", "jumpdist_backward", 1, reader_jumpdist }
};

static dofit_type_t	dofit_type = JUMPDIST_ALL;
static msl_fit_type_t	fit_type;

static void
usage(void)
{
	ecm_fprintf(ecm_stderr,
"Usage: lstatist [ <options> ] [<input file or dir>]\n"
" Options:\n"
"   -h: help\n"
"   -o: output file\n"
"   -t: type\n"
"   -T: model type\n"
	);
}

static void
parse_dofit_type(const char *typestr)
{
	ecm_uint32	i;

	for (i = 0; i < DOFIT_TYPE_MAX; i++) {
		if (ecm_strcmp(dofit_infos[i].desc, typestr) == 0) {
			dofit_type = (i + 1);
			return;
		}
	}

	ERR("%s: invalid target data", typestr);
	ecm_exit(3);
}

static void
parse_args(int argc, char *argv[])
{
	int	c;

	while ((c = ecm_getopt(argc, argv, "ho:t:T:")) != -1) {
		switch (c) {
		case 'h':
			usage();
			ecm_exit(0);
		case 'o':
			output = ecm_strdup(ecm_optarg);
			break;
		case 't':
			parse_dofit_type(ecm_optarg);
			break;
		case 'T':
			fit_type = parse_model_type(ecm_optarg);
			if (fit_type == MSL_FIT_NONE) {
				ERR("%s: invalid fit type", ecm_optarg);
				ecm_exit(2);
			}
			break;
		default:
			usage();
			ecm_exit(1);
		}
	}

	if (ecm_optind == argc) {
		usage();
		ecm_exit(1);
	}

	inputpath = ecm_strdup(argv[ecm_optind]);
}

static void
do_fitting(dofit_type_t type)
{
	msl_fit_t	*best_fit = NULL;
	ecm_double	goodness = -1;
	ecm_uint32	i;

	dofit_info_t	*pinfo = &dofit_infos[type - 1];
	char	*inputfile = NULL;
	
	if (ecm_is_directory(inputpath)) {
		char	*fname;

		ecm_asprintf(&fname, "log_%s.log", pinfo->logname);
		inputfile = ecm_gen_path(inputpath, fname, NULL);
		ecm_free(fname);
	}
	else
		inputfile = ecm_strdup(inputpath);

	if (!read_input(inputfile, pinfo->reader)) {
		ecm_free(inputfile);
		return;
	}
	ecm_free(inputfile);

	for (i = 0; i < MSL_FIT_MAX; i++) {
		msl_fit_t	*fit;
		ecm_double	goodness_this;

		if (fit_type != MSL_FIT_NONE && fit_type != (i + 1))
			continue;

		fit = lstat_fitting(pinfo->desc, i + 1, n_inputs, xvalues, yvalues);
		if (fit == NULL)
			continue;
		goodness_this = lstat_get_goodness(fit);
		if (goodness < 0 || goodness > goodness_this) {
			lstat_fitting_free(best_fit);
			goodness = goodness_this;
			best_fit = fit;
		}
	}

	if (best_fit == NULL) {
		ERR("fitting error");
	}	
	else {
		lstat_report_fitting(best_fit, pinfo->xdiff);
		printf("goodness: %lf\n", goodness);
	}
}

int
main(int argc, char *argv[])
{
	ECM_initialize();

	ecm_set_output_file(ecm_stderr);
	ecm_set_verbose_mode(ECM_TRUE);

	parse_args(argc, argv);

	do_fitting(dofit_type);

	return 0;
}
